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Abstract 

This paper reports the resuhs of an analysis of the Doppler tracking data of 
Pioneer probes which did show an anomalous behaviour. A software has been 
developed for the sake of performi ng a data analysis as i ndependent as possible 



from that of J. Anderson et al. (jAnderson et al.1 . 120021 ). using the same data 



set. A first output of this new analysis is a confirmation of the existence of 
a secular anomaly with an amplitude about 0.8 nms~^ compatible with that 
reported by Anderson et al. A second output is the study of periodic variations 
of the anomaly, which we characterize as functions of the azimuthal angle ip 
defined by the directions Sun-Earth Antenna and Sun-Pioneer. An improved fit 
is obtained with periodic variations written as the sum of a secular acceleration 
and two sinusoids of the angles ip and 2ip. The tests which have been performed 
for assessing the robustness of these results are presented. 

Key words: Pioneer anomaly, Doppler tracking data, interplanetary 
trajectography 



1. Introduction 

The Pioneer 10 and 11 spacecrafts performed accurate celestial mechanics 
experiment after their encounters with Jupiter and Saturn. Doppler tracking 
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by the NASA Deep Space Network ( Movei . 2000l ) of the two spacecrafts which 
followed escape orbits to opposite ends of the solar system revealed an anomalous 
behaviour which remains imperfectly understood. Precisely, the tracking data 
do not meet the expectations drawn from standard gravity force law and a better 



fit is obtained by adding a constant Sun-ward acceleration ([Anderson et al 
1998f) . This anomalous acceleration, now refered to as Pioneer anomaly, has a 



similar magnitude for the two spacecrafts. The valu e reported by the JPL team 
who discovered the anomaly is 0.874±0.133nms~^ ( Anderson et al.l . [2OO2I ). 

The exte nsive analysis by the J PL team has been published after years of 
cross checks ( Anderson et al.l 12002 ) . The presence of an ano maly and it s mag- 
nitude have been confirmed by different analysis software (jMarkwardtl . 12002 ; 
Olsenl I2OO7I ). A number of mechanisms have been considered as attempts of 



explanations of the anomaly as a systematic effect generated b y the spacecraft 
itself or its environment (see as an example ( Nieto et al. . 20051 )) but they have 
not led to a satisfactory understanding to date. If confirmed, the Pioneer sig- 
nal might reveal an anomalous behaviour of gravity at scales of the order of 
the size of the solar system and thus have a strong impact on fundamental 
physics, astrophysics and cosmology. It is therefore important to explore all its 
facets in the wider contex t of navigation and test of the gravity law i n the solar 
system (see for example ( Revnaud &: Jaekel . 20081 : Bertolami et al. . 20061 ) and 
references therein). 

Crucial i nformations are expected to come out from th e re-analysis of Pi- 
oneer data ( Nieto fc Anderson . 2005t Turvshev et al. . 2006f ). An international 
collaboration has been built t o this aim within the frame of International Space 
Science Institute (llSSj. l200d). The Pionee r data which had been analysed by 
the Anderson team ( Anderson et al. ■ l2002h have been made available by Slava 
Turyshev (JPL) in the framework of this collaboration. They consist in Orbit 
Data Files (ODF) which contain in particular Pioneer 10 Doppler data from 
November 30th, 1986 to July 20th, 1998. 

The aim of the present paper is to report some results of the analysis of 
these data performed by a collaboration between thr ee groups at ONERA, OCA 



and LKB within the "Groupe Anomalie Pioneer" (|GAPl . I2006h . A dedicated 
software called ODYSSEY has been developed to this purpose. The first result 
is to confirm the existence and magnitude of the anomalous secular acceleration 
reported by Anderson et al., using different and as independent as possible tools. 
The main motivation of the present paper is to study the periodic variations of 
the anomaly, which are known to exist besides the secular Pioneer anomaly. 

The presence of diurnal and seasonal variation s in the residuals has been dis- 
cussed by Anderson et al. ( Anderson et al. . 2002f) . We may emphasize that such 



modulated anomalies are unlikely to be due to anything relate d to the space- 
craft or its environment. Anderson et al. ( Anderson et al.l . l2002l) have proposed 
that they are due to modeling errors such as errors in the Earth's ephemeris, 
the orientation of the Earth's spin axis or the station's coordinates. However, 
these parameters are strongly constrained by other observational methods and 
it seems difficult to change them enough to explain the periodic anomaly. We 
will show here that this periodic anomaly can be at least partly represented in 
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terms of a modulation of the Doppler signal as a function of a unique azimuthal 
angle having a physical meaning. 



2. Development of the ODYSSEY software 

The Pioneer 10/11 spacecrafts were tracked by the Deep Space Network 
(DSN) antennas. A S-Band signal at about 2.11GHz is emitted at time ti by 
a DSN antenna and received onboard the spacecraft at time i2- The frequency 
is multiplied by a constant ratio of 240/221 by a transponder and sent back 
to a ground antenna where it was received at time t^. The Doppler shift, 
that is the difference between the up- and down- frequencies, is called a 2- 
way observable if the two antennas are the same, a 3-way observable if they 
are different. In fact, the ODF observable is the average of the Doppler shift 
over a time s pan called comp ression interval ([Mover . I2OOOI ). The ODF format, 
described in (jWacklevl . l2000f ). also contains the compression time, the date of 



the middle of the compression interval, the emitted frequency and the receiving 
and transmitting DSN antenna identifiers. 

To analyse the ODF data, a software called ODYSSEY has been devel- 
oped at the Observatoire de la Cote d'Azur within a collaboration with Onera. 
ODYSSEY stands for "Orbit Determination and phYsical Studies in the So- 
lar Environment Yonder". One of the motivation for the development of this 
software is the simulation of the expected effects of the Solar System Odyssey 
project. It is basically an interplanetary trajectory determination software. It 
performs numerical integration in rectangular coordinates of dynamical equa- 
tions to propagate the position and velocity of the spacecraft and variational 
equations to propagate the sensitivity of the position and velocity with respect 
to the initial conditions of position and velocity and other parameters to be 
fitted. The nu merical integration uses the Adams-Moulton-Cowell algorithm 
( Henricil . Il962h at an order 10. The use of this algorithm enables the direct 



integration of second order equations. 

The values of the parameters to be estimated are obtained through a best fit 
procedure. The difference between measured observables O and the calculated 
ones C, which depend on the parameters E to be fitted, is linearized in terms 
of variations of E. For N values of the observable, this corresponds to a linear 
system of N equations with as many variables as there are terms in E. This 
system is then solved in ODYSSEY with classical iterative least squares analysis. 

Maneuvers are taken into account as increments of velocity SV along the 
three directions. The dates of the maneuvers are provided by JPL and their 
amplitudes estimated as parameters in the best fit analysis. Constraints are 
imposed on the maximum values of the estimated SV in non radial directions. 

The dynamical model to compute the motion of the spacecraft includes grav- 
itational attraction by the main bodies of the solar system and direct radiation 
pressure. The motion of the spacecraft is computed using non relativistic grav- 
itational equations, which a re known to be sufhcient at the considered level of 



accuracy (IMarkwardtl l2002i) . For t he solar radiation pressure we use the same 



model as in ( Anderson et al. . 20021 ) with the reflectivity coefhcient kr = 1.77, 
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the mass m = 259kg of the spacecraft and a cross sectional area S = 5.9 
(orientation of the antenna supposed here to be constant). The accuracy of 
this model is not critical because the acceleration due to the solar radiation 
pressure has decreased to less than 0.2nms~^ for Pioneer 10 at the considered 
distances. Other known standard per turbations have been fo und to be negligible 



for Pioneer at the considered period (j Anderson et al.l . 120021 ). 



In its present status, the data analysis does not take into account the detailed 
thermal models of the spacecraft, currently under study by different groups. 
These models are expected to produce a slowly evoluting radiation force due to 
heat dissipation from the Radioisotope Thermoelectric Generators (RTG); this 
force should appear as a part of the anomalous secular acceleration to be found 
below. Improved analysis of this important issue will only be possible when 
accurate enough models of a time dependent radiative budget will be available. 

The dynamical equations are integrated in the Barycentric Celestial Refer- 
ence System (BCRS); The reference time scale is the Barycentric Coordinate 
Time (TCB). Positions of terrestrial stations are expressed in the International 
Terrestrial Reference System (ITRS); the reference time scale is the Coordi- 
nated Universal Time (UTC). The transition between ITRS and BCRS on the 
one hand, and TCB and UTC on the other hand, are performed according to 
2003 International Earth Rotation Service (lERS) conventions. The positions 
of celestial bodies are obtained from DE 405 ephemeris from JPL in BCRS with 
a reference time scale which is refered to as "Tep^" and which is similar to the 



Barycentric Dynamical Time ([Standishl . 119981 ) 



Special efforts have been devoted in the development of ODYSSEY for han- 
dling the calculation of the ODF observable. In a first step, the perturbations 
of the round-trip light time are not taken into account. The ODF observable, 
that is the average of the Doppler shift over the compression time, is computed 
through a numerical approximation using the 4-points Simpson method. For 
the available Pioneer data, the compression time does not exceed 1980s. The 
accuracy of the 4-points Simpson method in this case is evaluated to be 0.4 mHz. 

The instantaneous Doppler shift is calculated in te rms of velocities o f the 



endpoints, evaluated at the event times ti, t2 and ^3 (jMarkwardtl l2002i ). As 
only ^3 is provided in ODF, t2 and ti have to be determined, which is done 
iteratively using the light time equation 

^2=^3 0i23{t2,t3) , ti^t2 (3«12(tl,i2) (1) 

c c 

where S£i2 and 6i23 are respectively the Shapiro time delay affecting the up- 
and down-link signals. 

The spin frequency involved in Doppler shift estimation is provided by JPL. 
The ODF observable comes from JPL already corrected for data after the 17"^ 
of July 1990. Before this date, the spin correction is done by our software, 
taking the closest available value in the provided time series. As the maximal 
difference between two successive spin values in the file is 0.2 mHz, the accuracy 
of this method is considered as sufficient. 
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In the second step of the computation of the observable, the perturba- 
tions which affect the propagation of the tracking signal are taken into ac- 
count. The Shapiro ti me delay and the solar corona effect are modeled as in 
20021 ) . For the determination of the electron density necessary 



( Anderson et al 



to compute ionospheric effect, two models have been im plemented, the Inter- 
national Reference Ionosp here (IR I ) 200 7 (jBilitzal . 120011 ) and the Parametrized 
Ionospheric Model (PIM) (|Daniel]| . fl995h . 

Two models have been implemented in ODYSSEY for the m apping func tions 
of the tropospheric effect, the Niell Ma pping Functions (NM F) ( Niell . 19961 ) and 
the Global Mapping Function (GMF) (jBoehm et all l2006l ). 

All these perturbations are modeled as delays affecting the signal (except for 
the Shapiro delay, their effect in equation ([T]) is however negligible). They are 
written as added propagation lengths, d£i2 and (5^23 for the up- and down-links. 
Their effect on the ODF observable is therefore given by the difference of these 
quan tities between the "start" and "stop" of the compression interval (| Mover . 
20001) . 



3. Confirmation of the existence of a constant anomaly 

Our first aim was to study t he secular Pioneer anomaly reported by Anderson 



et al. ( Anderson et al.l . l2002h . To this aim, we performed a best-fit with a 



constant anomalous acceleration ap exerted on the probe. The acceleration was 
centered on a point chosen at the Sun, the Solar System Barycentre (SSB) or the 
Earth. The results obtained with these 3 possibilities were not distinguishable 
from each other, due to the large distance of the probe to the center of the solar 
system. From now on, we consider the center to be the Sun. 

The initial conditions as well as the three components of each maneuver 
are also fitted. The components of the maneuvers are fitted in the RTN (Ra- 
dial Transverse Normal) frame and the transverse and normal amplitudes are 
constrained to be inferior to 0.2 ms^^. The IRI 2007 model for ionospheric cor- 
rection and GMF for tropospheric correction are used. Points with an elevation 
inferior to 20° are rejected so as to limit the effect of imperfections of atmospher- 
ical models. Outliers are also rejected when their difference with the expectation 
exceeds 100 Hz at the first iteration and 6cr at the following iterations with cr 
the standard deviation of the residuals at this iteration. 

The analysis performed with the software ODYSSEY confirms that a better 
fit is obtained with a constant sunward acceleration. The value estimated by 
ODYSSEY for the anomalous acceleration is ap = 0.84 ± 0.01 nms~^ with the 
formal error given at la. The postfit residuals show a standard deviation of 
9.8 mHz, which is largely improved with respect to a fit without ap. These 
residuals are shown on Fig. [TJ 

On t he Fig. [21 is given anot her representation agreeing with figure 8 of An- 



derson (jAnderson et al.l 120021 ). This figure is reconstructed by the following 
procedure: the postfit residuals are expressed as O — C{Xq^* , Xq^* , M°p*' , ap) 
where Xq^* , Xq^*^ and M°p^ are the optimized values for initial position, ini- 
tial velocity and maneuver amplitudes; then the Fig. [2] represents the quantity 
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Figure 1: Best fit residuals of the Doppler 
tracking data of Pioneer 10 with an anoma- 
lous acceleration ap. 



Figure 2: Reconstruction of an anomalous 
contr ibution, as in Fig. 8 of llAnderson et all 
1 mms ^ corresponds to 15.4mHz. 



O — C (Xq^* , Xq^^ , M°P*', 0) where the anomalous acceleration has been nullified; 
this representation highlights the need of the constant acceleration to reduce 
the residuals. 

It can be emphasized that the level of the residuals on Fig. [T] is higher than 
the measurement noise. It is also clear on the figure that the postfit residuals do 
not correspond to a white gaussian noise. In order to highlight the existence of 
systematic structures in the residuals, we zoom on a time interval corresponding 
to the period from 23 November 1996 to 23 December 1996 where Pioneer 10 
was on opposition (Sun, Earth and Pioneer 10 aligned in this order). The data 
set is thus less affected by solar plasma and it shows daily variations of the 
residuals (see Fig. [3]). This figure clearly shows the existence of daily variations 
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Days from 23 November 1 992 



Figure 3: Best fit residuals of the Doppler tracking data of Pioneer 10, for a one-month period 
near opposition. Different symbols or colors refer to different couples of stations (sec the text 
for details). 

in the residuals, which are emphasized by distinguishing the data according to 
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different couples of emission/reception antennas. Different symbols or colors 
refer to different couples of stations, with M standing for Madrid (DSN antenna 
63), G for Goldstone (DSN antenna 14), C for Canberra (DSN antenna 43). The 
full line is a daily sinusoid fitted to the residuals of the dominant data (M,M). 



4. Study of the periodic variations of the Pioneer Anomaly 

In order to characterize the periodic variations of the residuals, we have per- 
formed a spectral analysis of the postfit residuals. As the Pioneer data points 
are not evenly distributed, the commonly used Fourier transform is not appro- 
priate for the analysis. We have use d the software Spar Spcc which estimates 



frequency distributions in time series (IBourguignon et al.l , ,2007 ). The result of 



this spectral analysis is shown on Fig. \E\ 

The presence of significant periodic terms is clear at the periods measured 
with respect to a day — 86400 s: 

/i = 0.9974 ±0.0004 day 
/2 ^ i(0.9972± 0.0004) day 
/3 = 189 ± 32 day. 

As 0.9972 day = 1.0 sidereal day, these periods are consistent with variations 
on one sidereal day, half a sidereal day, and half a year. 

The presence of diurnal an d seasonal variations in the residuals has also been 
reported by Anderson et al. ( Anderson et al.l . l2002f ). It has to be noted that 



these specific periods are unlikely to be due to anything related to the spacecraft 
or its environment. Errors in the atmospherical models would induce daily vari- 
ations. As these effects depend on the conditions at the stations, such errors 
would be expected at solar day rather than at sidereal day. Anderson et al. 
( Anderson et al. . 2002f l propose modeling errors such as in the case of Earth's 



ephemeris, the orientation of Earth's spin axis or the station's coordinates. How- 
ever, these parameters are strongly constrained by other observational methods 
and it seems difficult to change them enough to explain the periodic anomaly. 

The main motivation of the present paper is to test an alternative explanation 
where some perturbation would modify the propagation of the tracking signal 
along the path from the Earth antenna and the spacecraft. The idea is to 
represent such a perturbation, whatever its origin, as a function of the angle ip 
defined as the difference between the Earth Antenna (A) azimuthal angle and 
the Pioneer (P) azimuthal angle : (p = Lpp — (pA (see Fig.Hl). The main interest 
of this geometrical model is that it should simultaneously account for the orbital 
movement of the Earth around the Sun and the diurnal rotation of the Earth. 

As this perturbation is supposed to be periodic, it will be represented by a 
few Fourier coefficients 

A/ = ^ (u„(cos(n(p„) + cos{nipd)) + w'„(sin(ni^„) + sm{nipd))) (2) 
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Figure 4: Definition of the angle ip. 



Here ipu and >^d are the angles evaluated on the up- and down-links. Several 
tests have been performed with different numbers of coefficients Vn and v'^ fitted 
in addition to the constant anomalous acceleration model. The best results 
have been obtained with superpositions of functions with n = \ and 2. The 
addition of higher order coefficients has not modified significantly the residuals. 
To summarize, the new fit is identical to the previous one but for the addition 
of the following modification of the Doppler signal: 

A/ = wi(cos((y3„) -I- cos(v5d)) -I- v\[sv\i(ipu) + sin((^d)) 
+V2{cos{2ipu) + cos(2(^d)) -I- W2(sin(2(/3„) + sin(2(pd)) 

This model results in a spectacular improvement of the best fit residuals, 
with the standard deviation reduced by a factor of the order of two. Precisely, 
the standard deviation of the residuals which was 9.8 mHz without the peri- 
odic terms has been decreased to 5.5 mHz. The values of the fitted anomalous 
parameters are reported in table [T] with two different ionospheric models. 

The robustness of the results with respect to other modifications of the pro- 
cedure is evaluated by the few tests which follow. Using the IRI model and 
selecting data according to minimal elevation or outlier definition, we obtain 
the results shown in the tables [2] and [3l These tables show that the best fit 
results are robust enough, except for the amplitude v\ , when details of the pro- 
cedure are changed. In particular, the anomalous secular acceleration is not very 
sensitive to these changes and it remains roughly the same that in the previous 
best analysis without modulated terms. For the parameters other than v\, the 
uncertainty in the determination has to be estimated from the variations of the 
results when the details are changed, and certainly not from the formal errors 
given for a given fit which are smaller. 

For the amplitude ui, which is clearly less stable than the other parameters, 
it is better to consider that it is not determined by the best fit procedure. This 
feature can plausibly be understood from the following argument. The direction 
from the Sun to Pioneer f varies only slowly, because of the large distance of 
the probe, ft follows that the periodic terms contain essentially variations at 
one year and one sidereal day for n = 1 , at half a year and half a sidereal day for 
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Figure 5; SparSpec analysis of the residuals 
from the fit with a constant acceleration. 
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Figure 6: SparSpec analysis of the residuals 
from the fit with a constant acceleration and 
periodic terms. 



n = 2. The fit of the initial conditions also induces terms at the same periods. 
In particular, a change of the initial conditions may easily produce var iations 
masking modulated terms ( Courty et al. . 2008t Jaekel Sz Reynaudl - lioosh . This 
is probably the reason why the yearly period, which corresponds to a large 
potential amplitude, was in fact not detected in the spectral analysis of the best 
fit residuals with only a secular anomaly. 

Now the best fit with modulated terms tell us a different story. The modu- 
lated terms are looked for in a dedicated best fit procedure and they are unam- 
biguously found to differ from zero. The values of ap, v[, and v'2 seem to have 
robust estimates, while v\ is less stable. The analysis of the correlation shows 
a large value between v\ and the initial conditions (0.814 compared to 1 for a 
total correlation). There is an interdependance between the estimated values 
of v\ and the initial conditions. This large correlation reflects the phenomenon 
described in the previous paragraph explaining the instability of v\ estimation. 

An even more impressive demonstration of the improvement of the data 
analysis drawn by the inclusion of modulated terms comes from the spectral 
analysis of the residuals. This spectral analysis is represented on Fig. [6] for the 
best fit with modulated terms. It is drawn intentionnally at the same scale as 
Fig. [5] so that one can easily notice the global reduction of the main peaks in 
the spectrum as well as all the secondary peaks. 



5. Conclusion 



In the present paper, we have reported the first results of our re-analysis 
of the Pioneer 10 Doppler data for the 1 986 to 1998 time span. The data 
are the same that were already analysed in ( Anderson et al. , 2002 ; Markwardtl 
2002l : IOlsenl . l2nn7h but they have been dealt with by using a dedicated software 
ODYSSEY developed to this purpose. The improvement of the data fit with 
a constant anomalous acceleration exerted on Pioneer 10 has been confirmed 
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by this new data analysis. Its magnitude is compatible with that reported by 
Anderson et al. 

The paper has then been focused on the study of periodic terms in the residu- 
als, which were already noticed in previous studies and are clearly revealed using 
spectral analysis. The main new result of the paper is that a large part of these 
diurnal and seasonal anomalies may be captured in a simple geometrical model 
where the light time on the tracking path is modified in a manner depending 
only on the azimuthal angle Lp between the Sun-Earth and Sun-probe lines. This 
geometrical model could represent in a simple way the physical effects expected 
on light propagation in some metric extensions of general relativity which have 
been stud ied as potential candidates for the explanation of the secular Pioneer 
anomaly ( Jaekel fc Revnaud . 20051 20061 ). The search for modulated anomalies 
might even be used to distinguish those c andidates from other which affect much 
less the propagation of tracking signals (jMoffatl . l200fih . However, the results of 
the paper cannot be considered as pointing to a particular possible explanation 
of the anomaly. At the present stage of the data analysis, similar effects could 
for example be obtained through a mismodeling of the solar corona model. 

Anyway, considering the modulated anomalies has allowed us to reduce by a 
factor of the order of two the standard deviation of the residuals. The results of 
the best fit have been found to be robust enough with respect to changes in the 
tropospheric or ionospheric model, as well as variations of the details of the best 
fit procedure. The most impressive output of the new analysis is given by the 
SparSpec analysis of the best fit residuals which shows a reduction of all periodic 
structures which were present in the residuals of the best fit without modulated 
terms. It is worth emphasizing that this has been done by considering only 
a function of the geometrical angle ip, and not by using different explanations 
for the diurnal and seasonal anomalies. This suggests that the new analysis 
constitutes a richer characterization of the Pioneer data, now involving not only 
a secular acceleration but also modulated terms, which will have to be compared 
with any, existing as well as future, possible explanation of the anomaly. 
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Table 1: Results of the best fit with periodic terms in the signal, with two diflferent ionospheric 
models. 



TEC Model 


IRI 2007 


PIM 


ap (pms^'^) 


-836±1 


-836±1 


Vi (mHz) 


124.3±9.3 


141.6±9.3 


v'l (mHz) 


-125.3±0.6 


-127.3±0.6 


V2 (mHz) 


2.7±0.2 


3.0±0.2 


v'2 (mHz) 


-4.8±0.1 


-4.9±0.1 


Residuals (mHz) 


5.5 


5.5 



Tabic 2: Results of the best fit with periodic terms in the signal, when varying the minimal 



evaluation for the considered data. The initial number of points is 1980-5. 



Aliii. (dc\ati(jii 


0° 


10= 


20" 


30" 


Up (pms^"^) 


-844±1 


-843±1 


-836±1 


-826±1 


Vi (mHz) 


86.4±9.3 


98.0±9.2 


124±9.3 


95±10.5 


v{ (mHz) 


-111.8±0.6 


-119.5±0.7 


-125±0.6 


-127±0.6 


V2 (mHz) 


2.1±0.2 


2.3±0.2 


2.7±0.2 


2.1±0.2 


V2 (mHz) 


-4.7±0.1 


-4.7±0.1 


-4.8±0.1 


-4.6±0.1 


Nr of eliminated points 


789 


911 


2996 


8584 


Residuals (mHz) 


6.0 


5.9 


5.5 


5.1 



Table 3: Results of the best fit with periodic terms in the signal, when varying the outlier 
criterion at Ncr. The initial number of points is 19805. 



Outliers at Na 


N=3 


N=6 


N=8 


ap (pms^'^) 


-806±1 


-836±1 


-848±1 


Vi (mHz) 


197±7.2 


124±9.3 


1G2±10.0 


v[ (mHz) 


-126±0.5 


-125±0.6 


-124±0.7 


f2 (mHz) 


3.4±0.1 


2.7±0.2 


2.6±0.2 


v'2 (mHz) 


-5.0±0.09 


-4.7±0.1 


-4.7±0.1 


Nr of eliminated points 


3879 


2996 


2933 


Residuals (mHz) 


4.16 


5.5 


5.9 
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